Mengestimasi Integral dengan Monte Carlo Method

Metode Monte Carlo

Monte Carlo merupakan salah satu metode komputasi dengan mengandalkan pengambilan sampel acak secara berulang untuk mendapatkan hasil numerik yang mendekati nilai sebenarnya. Konsep yang mendasari dari metode ini dengan menggunakan nilai acak untuk memecahkan masalah deterministik dalam prinsipnya. Metode Monte Carlo sering digunakan dalam memecahkan masalah matematika dan fisika dan peling berguna ketika sulit atau tidak mungkin untuk menggunakan pendekatan lain.

Dalam perkiraan integral dengan dengan motode Monte Carlo untuk memperkirakan area dan volume. Kita gunakan $$ \int_{0}^1 f(x)dx \approx {1 \over n} \sum_{i=0}^{n} f(x_i) $$

\int_{0}^{1} \int_{0}^{1} \int_{0}^{1} f(x,y,z)dx \space dy \space dz \approx {1 \over n} \sum_{i=0}^{n} f(x_i, y_i, z_i)

dimana $ x_i $ merupakan urutan angka acak dalam rentang satuan dan $ (x_i, y_i, z_i) $ merupakan urutan acak unit N point cube. Namun pada umumnya $ \int_{0}^{1} f \approx$ (ukuran A) $ \times $ (rata-rata f lebih dari n poin acak dalam A).

Tugas Programming

  1. Write program to verify numerically that $ \pi = \int_{0}^{2} (4 - x^2)^{1/2} $. Use the Monte Carlo method and 2500 random numbers.
  2. Use the Monte Carlo method to approximate the integral. $ \int_{-1}^{1} \int_{-1}^{1} \int_{-1}^{1} (x^2 + y^2 + z^2) dx \space dy \space dz $

Berikut merupakan program untuk mengestimasi integral menggunakan Monte Carlo method dengan bahasa pemrograman python.

Program berikut untuk melakukan pengujian bahwa $ \pi = \int_{0}^{2} (4 - x^2)^{1/2} $, dengan angka acak sebanyak 2500.
import random as r

analitik = 3.14159

count = 2500
x0 = 0; x1 = 2

def monte_carlo(f, count, x0, x1):
    in_area = 0.0
    for i in range(count):
        x_ = r.uniform(x0, x1);   
        if x0 <= x_ <= x1:
            in_area += f(x_)

    area_box = (x1 - x0)
    return ((in_area / count) * area_box)

def f(x):
    return ((4 - x**2)**(1/2))

hasil = monte_carlo(f, count, x0, x1)
print('N =', count)
print('Hasil aproksimasi:', hasil)
print('Relative true error:', (analitik - hasil)/analitik, '|', str(round(((analitik - hasil)/analitik)*100, 2)) + '%')

Dalam program tersebut, terdapat library random untuk mendapatkan nilai acak. program tersebut terdapat variabel analitik yang berisi nilai perhitungan secara analitik/nilai sebernarnya. Kemudian terdapat variabel count sebagai N (banyaknya perulangan), variabel x0 untuk batas bawah dan x1 untuk batas atas integral. Didalam fungsi monte_carlo(param1,param2..) tersebut, terdapat paramater f untuk fungsi integral yang akan digunakan, count, x0, dan x1. Kemudian dalam fungsi tersebut terdapat variable in_area, dimana variabel tersebut akan digunakan untuk menampung nilai acak yang dihasilkan pada area dari persamaan tersebut. kemudian terdapat looping yang akan mengulang sebanyak N. Dalam perulangan tersebut, tedapat variabel x_ untuk nilai acak dalam rentang x0 dan x1. Kemudan akan akan dilakukan pengecekan kembali apabila nilai dari x_ dalam rentang x0 dan x1, maka akan menjumlah terus menerus nilai dari perhitungan f(x), dimana ilustrasinya nilai tersebut merupakan titik dalam bentang persamaan yang didefinisikan. Setelah melakukan looping sebanyak N, maka akan mencari nilai dari area_box, area_box merupakan area dalam interval x0 dan x1 yang berbentuk kotak. Kemudian fungsi tersebut akan mengembalikan nilai dari perhitungan rata-rata dari in_area dikali dengan area_box.

Berikut merupakan ilustrasi dari program tersebut dengan persamaan $ f(x) = (4 - x^2)^{1/2} $.

.

Berikut merupakan hasil running program, dalam running program tersebut digunakan N sebanyak 2500 untuk running pertama dengan hasil perkiraan, yaitu 3.1169340168166424 dengan relative true error sebesar 0,78%. Kemudian dalam running kedua digunakan N sebanya 5000 dengan hasil perkiraan 3.1415150955638516. Nilai perkiraan yang dihasilkan akan dinamis setiap program dijalankan, terkadang error lebih besar atau lebih kecil. Namun, dengan memberikan nilai N yang lebih besar error yang dihasilkan dominan lebih kecil.

N = 2500
Hasil aproksimasi: 3.1169340168166424
Relative true error: 0.007848249829977014 | 0.78%

N = 5000
Hasil aproksimasi: 3.1415150955638516
Relative true error: 2.3842842684222457e-05 | 0.0%
Selanjutnya program kedua untuk mengaproksimasi nilai integral $ \int_{-1}^{1} \int_{-1}^{1} \int_{-1}^{1} (x^2 + y^2 + z^2) dx \space dy \space dz $
import random as r

analitik = 8.0

count = 1000
x0 = -1; x1 = 1; y0 = -1; y1 = 1; z0 = -1; z1 = 1

def monte_carlo(f, count, x0, x1, y0, y1, z0, z1):
    in_area = 0.0
    for i in range(count):
        x_ = r.uniform(x0, x1); y_ = r.uniform(x0, x1); z_ = r.uniform(z0, z1);
        if (x0 <= x_ <= x1) and (y0 <= y_ <= y1) and (z0 <= z_ <= z1):
            in_area += f(x_, y_, z_)

    area = (x1 - x0) * (y1 - y0) * (z1 - z0)
    return (in_area / count) * area

def f(x, y, z):
    return (x**2 + y**2 + z**2)

hasil = monte_carlo(f, count, x0, x1, y0, y1, z0, z1)
print('N =', count)
print('Hasil aproksimasi:', hasil)
print('Relative true error:', (analitik - hasil)/analitik, '|', str(round(((analitik - hasil)/analitik)*100, 2)) + '%')

Dalam program tersebut hampir sama dengan program pertama, hanya saja dalam program tersebut untuk mengaproksimasi integral lapis 3. Kemudian x0 dan x1 untuk integral pertama, y0 dan y1 untuk integral kedua, dan z0 dan z1 untuk integral ketiga. Proses perhitungan dan penyelesaian yang dilakukan hampir sama dengan program pertama, yaitu dengan menggunakan angka random dalam rentang yang ditentukan, dan titiknya juga berada pada dimensi yang berbeda.

Berikut merupakan hasil running program, dalam running program pertama digunakan N sebanyak 1000, maka dihasilkan nilai aproksimasi, yaitu 7.918782014912802 dengan relative true error sebesar 1.02%. Kemudian dalam running yang kedua menggunakan N sebanyak 5000, maka hasil aproksimasi yaitu 7.974525782192673 dengan relative true error sebesar 0.32%. Nilai perkiraan yang dihasilkan akan dinamis setiap program dijalankan, terkadang error lebih besar atau lebih kecil. Namun, dengan memberikan nilai N yang lebih besar error yang dihasilkan dominan lebih kecil.

N = 1000
Hasil aproksimasi: 7.918782014912802
Relative true error: 0.010152248135899766 | 1.02%

N = 5000
Hasil aproksimasi: 7.974525782192673
Relative true error: 0.0031842772259158547 | 0.32%
Berikut merupakan proses analitiks dari integral tersebut, dimana hasil analitiksnya yaitu 8. $$ I = \int_{-1}^{1} \int_{-1}^{1} \int_{-1}^{1} (x^2 + y^2 + z^3) dx \space dy \space dz \\ I = \int_{-1}^{1} \int_{-1}^{1} \left[ {1 \over 3}x^3 + y^2 + z^2 \right]_{-1}^{1} dy \space dz \\ I = \int_{-1}^{1} \int_{-1}^{1} \left[ {1 \over 3}(1)x^3 + (1)y^2 + (1)z^2 \right] - \left[ {1 \over 3}(-1)x^3 + (-1)y^2 + (-1)z^2 \right] \\ I = \int_{-1}^{1} \int_{-1}^{1} \left({2 \over 3} + 2y^2 + 2z^2 \right) dy \space dz \\ I = \int_{-1}^{1} \left[ {2 \over 3}y + {2 \over 3}y^3 + 2z^2y \right]_{-1}^{1} dz \\ I = \int_{-1}^{1} \left[ {2 \over 3}(1) + {2 \over 3}(1)^3 + 2z^2(1) \right] - \left[ {2 \over 3}(-1) + {2 \over 3}(-1)^3 + 2z^2(-1) \right] \\ I = \int_{-1}^{1} \left({8 \over 3} + 4z \right) dz \\ I = \left[{8 \over 3}z + {4 \over 3}z^3 \right]_{-1}^{1} \\ =\left[ {8 \over 3}(1) + {4 \over 3}(1)^3 \right] - \left[ {8 \over 3}(-1) + {4 \over 3}(-1)^3 \right] \\ = {24 \over 3} \\ = 8 $$

Sekian Terimakasih.